from scipy import eye
from scipy.integrate import ode
from scipy import zeros

y0, t0 = [1.0, 2.0, 3.0], 0
arg = 2;

def f(t, y):
    global arg
    dy = zeros(3,float)
    dy[0] = arg*y[0] + y[1]
    dy[1] = -arg*y[1]**2
    dy[2] = arg*y[2]/2-y[1]
    return [dy[0], dy[1], dy[2] ]
r = ode(f).set_integrator('zvode', method='bdf', with_jacobian=False, rtol=0.1,atol = 0.5)
r.set_initial_value(y0, t0)
t1 = 10
dt = 1
while r.successful() and r.t < t1:
    r.integrate(r.t+dt)
    print r.t, r.y

